singular value decompositions \(\mathbf{A} = \mathbf{U} \Sigma \mathbf{V}^T\)
iterative methods for numerical linear algebra
Except for the iterative methods, most of these numerical linear algebra tasks are implemented in the BLAS and LAPACK libraries. They form the building blocks of most statistical computing tasks (optimization, MCMC).
Our major goal (or learning objectives) is to
know the complexity (flop count) of each task
be familiar with the BLAS and LAPACK functions (what they do)
do not re-invent wheels by implementing these dense linear algebra subroutines by yourself
understand the need for iterative methods
apply appropriate numerical algebra tools to various statistical problems
All high-level languages (Julia, Matlab, Python, R) call BLAS and LAPACK for numerical linear algebra.
Julia offers more flexibility by exposing interfaces to many BLAS/LAPACK subroutines directly. See documentation.
2 BLAS
BLAS stands for basic linear algebra subprograms.
See netlib for a complete list of standardized BLAS functions.
Matlab uses Intel’s MKL (mathematical kernel libaries). MKL implementation is the gold standard on market. It is not open source but the compiled library is free for Linux and MacOS. However, not surprisingly, it only works on Intel CPUs.
Julia uses OpenBLAS. OpenBLAS is the best cross-platform, open source implementation. With the MKL.jl package, it’s also very easy to use MKL in Julia.
Typical BLAS functions support single precision (S), double precision (D), complex (C), and double complex (Z).
3 Examples
The form of a mathematical expression and the way the expression should be evaluated in actual practice may be quite different.
Some operations appear as level-3 but indeed are level-2.
Example 1. A common operation in statistics is column scaling or row scaling \[
\begin{eqnarray*}
\mathbf{A} &=& \mathbf{A} \mathbf{D} \quad \text{(column scaling)} \\
\mathbf{A} &=& \mathbf{D} \mathbf{A} \quad \text{(row scaling)},
\end{eqnarray*}
\] where \(\mathbf{D}\) is diagonal. For example, in generalized linear models (GLMs), the Fisher information matrix takes the form \[
\mathbf{X}^T \mathbf{W} \mathbf{X},
\] where \(\mathbf{W}\) is a diagonal matrix with observation weights on diagonal.
Column and row scalings are essentially level-2 operations!
usingBenchmarkTools, LinearAlgebra, RandomRandom.seed!(257) # seedn =2000A =rand(n, n) # n-by-n matrixd =rand(n) # n vectorD =Diagonal(d) # diagonal matrix with d as diagonal
# This works only when Matlab is installedusingMATLABmat"""d = rand(5, 1)diag(d)"""
ArgumentError: ArgumentError: Package MATLAB not found in current path.
- Run `import Pkg; Pkg.add("MATLAB")` to install the MATLAB package.
Example 2. Innter product between two matrices \(\mathbf{A}, \mathbf{B} \in \mathbb{R}^{m \times n}\) is often written as \[
\text{trace}(\mathbf{A}^T \mathbf{B}), \text{trace}(\mathbf{B} \mathbf{A}^T), \text{trace}(\mathbf{A} \mathbf{B}^T), \text{ or } \text{trace}(\mathbf{B}^T \mathbf{A}).
\] They appear as level-3 operation (matrix multiplication with \(O(m^2n)\) or \(O(mn^2)\) flops).
Random.seed!(123)n =2000A, B =randn(n, n), randn(n, n)# slow way to evaluate tr(A'B): 2mn^2 flops@benchmarktr(transpose($A) *$B)
BenchmarkTools.Trial: 62 samples with 1 evaluation.
Range (min … max): 50.385 ms … 147.322 ms┊ GC (min … max): 0.00% … 0.00%
Time (median): 66.011 ms ┊ GC (median): 0.00%
Time (mean ± σ): 81.664 ms ± 31.991 ms┊ GC (mean ± σ): 0.58% ± 1.40%
█ ▁
█▅█▁▅▇▁▁▁▅▁▅▇▁▅▁▁▅▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▅▁▁▁▅▅▅▁▅▅▁▁▇▅▇▇▇▇▇▇▁▇▁▇▁▅ ▁
50.4 msHistogram: log(frequency) by time 130 ms <
Memory estimate: 30.52 MiB, allocs estimate: 2.
But \(\text{trace}(\mathbf{A}^T \mathbf{B}) = <\text{vec}(\mathbf{A}), \text{vec}(\mathbf{B})>\). The latter is level-1 BLAS operation with \(O(mn)\) flops.
# smarter way to evaluate tr(A'B): 2mn flops@benchmarkdot($A, $B)
BenchmarkTools.Trial: 2698 samples with 1 evaluation.
Range (min … max): 1.793 ms … 2.206 ms┊ GC (min … max): 0.00% … 0.00%
Time (median): 1.836 ms ┊ GC (median): 0.00%
Time (mean ± σ): 1.850 ms ± 48.023 μs┊ GC (mean ± σ): 0.00% ± 0.00%
▂▅███▇▇▇▆▆▆▅▅▅▃▃▂▂ ▂
▃▃▆█████████████████████▇▆▆▆▇▅▅▆▅▆▅▅▅▇▅▅▆▆▅▆▅▆▆█▆▆▆▅▆▇▅▆▄▆ █
1.79 msHistogram: log(frequency) by time 2.05 ms <
Memory estimate: 0 bytes, allocs estimate: 0.
Example 3. Similarly \(\text{diag}(\mathbf{A}^T \mathbf{B})\) can be calculated in \(O(mn)\) flops.
# slow way to evaluate diag(A'B): O(n^3)@benchmarkdiag(transpose($A) *$B)
BenchmarkTools.Trial: 93 samples with 1 evaluation.
Range (min … max): 50.588 ms … 75.902 ms┊ GC (min … max): 0.00% … 0.00%
Time (median): 52.324 ms ┊ GC (median): 0.00%
Time (mean ± σ): 53.967 ms ± 4.768 ms┊ GC (mean ± σ): 0.89% ± 1.54%
█ ▁ ▂
█▇█▇▆█▃▃▄▃▃▁▁▃▃▃▃▄▃▃▃▁▃▃▁▁▁▁▁▃▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃ ▁
50.6 ms Histogram: frequency by time 74 ms <
Memory estimate: 30.53 MiB, allocs estimate: 3.
# smarter way to evaluate diag(A'B): O(n^2)@benchmarkDiagonal(vec(sum($A .*$B, dims =1)))
BenchmarkTools.Trial: 1257 samples with 1 evaluation.
Range (min … max): 3.330 ms … 15.680 ms┊ GC (min … max): 0.00% … 9.60%
Time (median): 3.640 ms ┊ GC (median): 0.00%
Time (mean ± σ): 3.972 ms ± 774.977 μs┊ GC (mean ± σ): 9.38% ± 12.67%
▄█▅▃▁
▅███████▇▇▅▅▄▃▃▃▂▂▂▂▂▂▁▂▂▁▁▂▁▂▃▅▆▃▄▄▄▄▄▃▃▂▃▃▂▂▂▂▂▂▂▂▂▃▃▃▃▂▂ ▃
3.33 ms Histogram: frequency by time 5.93 ms <
Memory estimate: 30.53 MiB, allocs estimate: 5.
To get rid of allocation of intermediate arrays at all, we can just write a double loop or use dot function.
functiondiag_matmul!(d, A, B) m, n =size(A)@assertsize(B) == (m, n) "A and B should have same size"fill!(d, 0)for j in1:n, i in1:m d[j] += A[i, j] * B[i, j]endDiagonal(d)endd =zeros(eltype(A), size(A, 2))@benchmarkdiag_matmul!($d, $A, $B)
BenchmarkTools.Trial: 1465 samples with 1 evaluation.
Range (min … max): 3.319 ms … 13.782 ms┊ GC (min … max): 0.00% … 0.00%
Time (median): 3.374 ms ┊ GC (median): 0.00%
Time (mean ± σ): 3.409 ms ± 509.119 μs┊ GC (mean ± σ): 0.00% ± 0.00%
▅▄ ▁▃▂▃▇▇▇▇██▆▃▆▄▄ ▁▂
██▇███████████████▇██▆▆▆▄▅▅▅▅▅▄▃▃▃▃▃▃▂▃▃▂▁▂▂▁▁▂▁▁▁▁▂▂▁▁▂▁▂▂ ▄
3.32 ms Histogram: frequency by time 3.58 ms <
Memory estimate: 0 bytes, allocs estimate: 0.
Exercise: Try @turbo (SIMD) and @tturbo (SIMD) from LoopVectorization.jl package.
4 Memory hierarchy and level-3 fraction
Key to high performance is effective use of memory hierarchy. True on all architectures.
Flop count is not the sole determinant of algorithm efficiency. Another important factor is data movement through the memory hierarchy.
In Julia, we can query the CPU topology by the Hwloc.jl package. For example, this laptop runs an Apple M2 Max chip with 4 efficiency cores and 8 performance cores.
For example, Xeon X5650 CPU has a theoretical throughput of 128 DP GFLOPS but a max memory bandwidth of 32GB/s.
Can we keep CPU cores busy with enough deliveries of matrix data and ship the results to memory fast enough to avoid backlog?
Answer: use high-level BLAS as much as possible.
A distinction between LAPACK and LINPACK (older version of R uses LINPACK) is that LAPACK makes use of higher level BLAS as much as possible (usually by smart partitioning) to increase the so-called level-3 fraction.
To appreciate the efforts in an optimized BLAS implementation such as OpenBLAS (evolved from GotoBLAS), see the Quora question, especially the video. Bottomline is
Get familiar with (good implementations of) BLAS/LAPACK and use them as much as possible.
5 Effect of data layout
Data layout in memory affects algorithmic efficiency too. It is much faster to move chunks of data in memory than retrieving/writing scattered data.
Storage mode: column-major (Fortran, Matlab, R, Julia) vs row-major (C/C++).
Cache line is the minimum amount of cache which can be loaded and stored to memory.
x86 CPUs: 64 bytes
ARM CPUs: 32 bytes
In Julia, we can query the cache line size by Hwloc.jl.
# Apple Silicon (M1/M2 chips) don't have L3 cacheHwloc.cachelinesize()
ErrorException: Your system doesn't seem to have an L3 cache.
Accessing column-major stored matrix by rows (\(ij\) looping) causes lots of cache misses.
Take matrix multiplication as an example \[
\mathbf{C} \gets \mathbf{C} + \mathbf{A} \mathbf{B}, \quad \mathbf{A} \in \mathbb{R}^{m \times p}, \mathbf{B} \in \mathbb{R}^{p \times n}, \mathbf{C} \in \mathbb{R}^{m \times n}.
\] Assume the storage is column-major, such as in Julia. There are 6 variants of the algorithms according to the order in the triple loops.
jki or kji looping:
# inner most loopfor i in1:m C[i, j] = C[i, j] + A[i, k] * B[k, j]end
- `ikj` or `kij` looping:
# inner most loop for j in1:n C[i, j] = C[i, j] + A[i, k] * B[k, j]end
ijk or jik looping:
# inner most loop for k in1:p C[i, j] = C[i, j] + A[i, k] * B[k, j]end
We pay attention to the innermost loop, where the vector calculation occurs. The associated stride when accessing the three matrices in memory (assuming column-major storage) is
Variant
A Stride
B Stride
C Stride
\(jki\) or \(kji\)
Unit
0
Unit
\(ikj\) or \(kij\)
0
Non-Unit
Non-Unit
\(ijk\) or \(jik\)
Non-Unit
Unit
0
Apparently the variants \(jki\) or \(kji\) are preferred.
""" matmul_by_loop!(A, B, C, order)Overwrite `C` by `A * B`. `order` indicates the looping order for triple loop."""functionmatmul_by_loop!(A::Matrix, B::Matrix, C::Matrix, order::String) m =size(A, 1) p =size(A, 2) n =size(B, 2)fill!(C, 0)if order =="jki"@inboundsfor j =1:n, k =1:p, i =1:m C[i, j] += A[i, k] * B[k, j]endendif order =="kji"@inboundsfor k =1:p, j =1:n, i =1:m C[i, j] += A[i, k] * B[k, j]endendif order =="ikj"@inboundsfor i =1:m, k =1:p, j =1:n C[i, j] += A[i, k] * B[k, j]endendif order =="kij"@inboundsfor k =1:p, i =1:m, j =1:n C[i, j] += A[i, k] * B[k, j]endendif order =="ijk"@inboundsfor i =1:m, j =1:n, k =1:p C[i, j] += A[i, k] * B[k, j]endendif order =="jik"@inboundsfor j =1:n, i =1:m, k =1:p C[i, j] += A[i, k] * B[k, j]endendendusingRandomRandom.seed!(123)m, p, n =2000, 100, 2000A =rand(m, p)B =rand(p, n)C =zeros(m, n);
BenchmarkTools.Trial: 86 samples with 1 evaluation.
Range (min … max): 57.729 ms … 70.826 ms┊ GC (min … max): 0.00% … 0.00%
Time (median): 58.222 ms ┊ GC (median): 0.00%
Time (mean ± σ): 58.433 ms ± 1.412 ms┊ GC (mean ± σ): 0.00% ± 0.00%
▁ █ ▆ ▁▁▁▁ ▁ ▃ ▁ ▆ ▁
▇█▇▁█▇▄█▇▇████▄▇▁█▇▄▇█▄▇▄▄▁▇█▁▁█▄▄▁▁█▄▄▄▄▁▇▁▁▁▁▁▄▄▁▄▁▁▁▇▄▁▄ ▁
57.7 ms Histogram: frequency by time 59.3 ms <
Memory estimate: 0 bytes, allocs estimate: 0.
@benchmarkmatmul_by_loop!($A, $B, $C, "kji")
BenchmarkTools.Trial: 27 samples with 1 evaluation.
Range (min … max): 183.530 ms … 212.516 ms┊ GC (min … max): 0.00% … 0.00%
Time (median): 184.559 ms ┊ GC (median): 0.00%
Time (mean ± σ): 186.442 ms ± 5.699 ms┊ GC (mean ± σ): 0.00% ± 0.00%
█▆▄
███▁▁▆▄▆▁▁▁▄▁▁▁▄▁▁▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄ ▁
184 ms Histogram: frequency by time 213 ms <
Memory estimate: 0 bytes, allocs estimate: 0.
\(ikj\) and \(kij\) looping:
@benchmarkmatmul_by_loop!($A, $B, $C, "ikj")
BenchmarkTools.Trial: 10 samples with 1 evaluation.
Range (min … max): 509.252 ms … 527.728 ms┊ GC (min … max): 0.00% … 0.00%
Time (median): 515.239 ms ┊ GC (median): 0.00%
Time (mean ± σ): 515.454 ms ± 5.554 ms┊ GC (mean ± σ): 0.00% ± 0.00%
█ ▁ ▁ ▁ ▁ ▁ ▁▁ ▁
█▁▁▁▁▁█▁▁▁█▁▁▁▁▁█▁▁▁▁▁█▁▁█▁▁▁▁██▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
509 ms Histogram: frequency by time 528 ms <
Memory estimate: 0 bytes, allocs estimate: 0.
@benchmarkmatmul_by_loop!($A, $B, $C, "kij")
BenchmarkTools.Trial: 10 samples with 1 evaluation.
Range (min … max): 507.229 ms … 530.723 ms┊ GC (min … max): 0.00% … 0.00%
Time (median): 511.375 ms ┊ GC (median): 0.00%
Time (mean ± σ): 513.003 ms ± 6.944 ms┊ GC (mean ± σ): 0.00% ± 0.00%
████ ██ █ ██ █
████▁▁▁▁▁█▁█▁▁▁▁▁█▁▁██▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
507 ms Histogram: frequency by time 531 ms <
Memory estimate: 0 bytes, allocs estimate: 0.
\(ijk\) and \(jik\) looping:
@benchmarkmatmul_by_loop!($A, $B, $C, "ijk")
BenchmarkTools.Trial: 21 samples with 1 evaluation.
Range (min … max): 244.667 ms … 265.676 ms┊ GC (min … max): 0.00% … 0.00%
Time (median): 247.189 ms ┊ GC (median): 0.00%
Time (mean ± σ): 249.187 ms ± 5.163 ms┊ GC (mean ± σ): 0.00% ± 0.00%
█▃▃ ▃ ▃
▇███▁▁█▇▇█▁▁▁▁▁▇▁▁▁▁▁▁▇▁▇▁▇▁▇▁▁▁▁▇▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▇ ▁
245 ms Histogram: frequency by time 266 ms <
Memory estimate: 0 bytes, allocs estimate: 0.
@benchmarkmatmul_by_loop!($A, $B, $C, "ijk")
BenchmarkTools.Trial: 21 samples with 1 evaluation.
Range (min … max): 245.461 ms … 255.071 ms┊ GC (min … max): 0.00% … 0.00%
Time (median): 247.440 ms ┊ GC (median): 0.00%
Time (mean ± σ): 248.889 ms ± 3.108 ms┊ GC (mean ± σ): 0.00% ± 0.00%
█▁▁ ▁▁▁▁█▁▁ ▁ ▁ ▁ ▁ ▁ ▁ ▁▁ ▁
███▁▁▁█████▁██▁▁▁▁▁▁▁▁▁█▁▁█▁▁▁▁█▁▁▁▁▁█▁▁▁▁▁▁▁█▁█▁▁██▁▁▁▁▁▁▁▁█ ▁
245 ms Histogram: frequency by time 255 ms <
Memory estimate: 0 bytes, allocs estimate: 0.
Question: Can our loop beat BLAS? Julia wraps BLAS library for matrix multiplication. We see BLAS library wins hands down (multi-threading, Strassen algorithm, higher level-3 fraction by block outer product).
@benchmarkmul!($C, $A, $B)
BenchmarkTools.Trial: 1708 samples with 1 evaluation.
Range (min … max): 2.526 ms … 24.014 ms┊ GC (min … max): 0.00% … 0.00%
Time (median): 2.643 ms ┊ GC (median): 0.00%
Time (mean ± σ): 2.926 ms ± 1.033 ms┊ GC (mean ± σ): 0.00% ± 0.00%
██▆▂▃▂▂▃▃▁▁▁ ▁ ▁
████████████████▇▇█▇▇▇▆▆▅▆▄▄▆▆▁▁▆▄▄▅▁▁▄▄▄▁▁▄▄▁▁▁▄▁▅▁▁▁▄▄▄▄ █
2.53 msHistogram: log(frequency) by time 6.75 ms <
Memory estimate: 0 bytes, allocs estimate: 0.
# direct call of BLAS wrapper function@benchmarkLinearAlgebra.BLAS.gemm!('N', 'N', 1.0, $A, $B, 0.0, $C)
BenchmarkTools.Trial: 1653 samples with 1 evaluation.
Range (min … max): 2.554 ms … 26.166 ms┊ GC (min … max): 0.00% … 0.00%
Time (median): 2.672 ms ┊ GC (median): 0.00%
Time (mean ± σ): 3.022 ms ± 1.278 ms┊ GC (mean ± σ): 0.00% ± 0.00%
██▄▂▂▃▂▂▂▂▁ ▁
████████████████▅▇▆▆▆▆▅▇▆▇▄▄▄▅▄▄▁▄▄▄▄▄▁▄▄▄▄▄▁▄▁▄▁▄▄▄▁▄▄▄▁▄ █
2.55 msHistogram: log(frequency) by time 7.47 ms <
Memory estimate: 0 bytes, allocs estimate: 0.
Question (again): Can our loop beat BLAS?
Exercise: Annotate the loop in matmul_by_loop! by @turbo and @tturbo (multi-threading) and benchmark again.
6 BLAS in R
Tip for R users. Standard R distribution from CRAN uses a very out-dated BLAS/LAPACK library.
usingRCallR"""sessionInfo()"""
RObject{VecSxp}
R version 4.3.2 (2023-10-31)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Sonoma 14.4.1
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
locale:
[1] C
time zone: America/Los_Angeles
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
loaded via a namespace (and not attached):
[1] compiler_4.3.2
┌ Warning: RCall.jl:
│ Attaching package: 'dplyr'
│
│ The following objects are masked from 'package:stats':
│
│ filter, lag
│
│ The following objects are masked from 'package:base':
│
│ intersect, setdiff, setequal, union
│
└ @ RCall /Users/huazhou/.julia/packages/RCall/dDAVd/src/io.jl:172
┌ Warning: RCall.jl: Warning: Some expressions had a GC in every iteration; so filtering is disabled.
└ @ RCall /Users/huazhou/.julia/packages/RCall/dDAVd/src/io.jl:172
# A tibble: 1 x 13
expression min median `itr/sec` mem_alloc `gc/sec` n_itr n_gc
<bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl> <int> <dbl>
1 A %*% B 124ms 125ms 7.75 30.5MB 7.75 4 4
total_time result memory time
<bch:tm> <list> <list> <list>
1 516ms <dbl [2,000 x 2,000]> <Rprofmem [1 x 3]> <bench_tm [4]>
gc
<list>
1 <tibble [4 x 3]>
Re-build R from source using OpenBLAS or MKL will immediately boost linear algebra performance in R. Google build R using MKL to get started. Similarly we can build Julia using MKL.
Matlab uses MKL. Usually it’s very hard to beat Matlab in terms of linear algebra.
usingMATLABmat"""f = @() $A * $B;timeit(f)"""
ArgumentError: ArgumentError: Package MATLAB not found in current path.
- Run `import Pkg; Pkg.add("MATLAB")` to install the MATLAB package.
7 Avoid memory allocation: some examples
7.1 Transposing matrix is an expensive memory operation
In R, the command
t(A) %*% x
will first transpose A then perform matrix multiplication, causing unnecessary memory allocation